Journal of Computational Biology

The study of evolution between species (phylogenetics) and the study of evolution within a species (population genetics) are highly related, as the same biological mechanisms are fundamental to both fields. Although both have been studied for a long time, their joint study in a unified setting has been prevented by the different time scales they consider and the different data types they employ. A similar discrepancy holds for their whole-genome specializations, comparative genomics and pangenomics. Two active areas in these fields are genome rearrangement studies and graphical pangenomics, respectively. Since the emergence of graphical pangenomics, these have existed as separate fields, despite observations that central data structures representing genomic variants in both fields are highly similar. While there exists a wealth of theoretical results for various rearrangement models in comparative genomics, the application to pangenomic data is hampered by the limitations of rearrangement problem formulations. On the practical side, pangenomes typically contain too many individual genomes for classical problems, such as the often NP-hard parsimony problems, to be solved, or for all-vs-all comparisons using rearrangement distances to be performed. On the theoretical side, some assumptions in the formulation of rearrangement problems, such as the assumption of an underlying tree, are inadequate for many pangenomes. In this work, we propose the Complete Ancestral Reconstruction for Pangenomes (CARP) problem, which overcomes these limitations while retaining intuitive relationships to both classical rearrangement problems and pangenome graphs.

Integrative analysis of the properties of multiple networks pertaining to a biological system is a key problem in systems biology. But computing the properties of thousands of large networks poses a challenge. Most techniques that partially address this challenge, such as parallelism and optimized algorithms, treat the analysis of each network as an independent task. However, the networks arising from the same biological system are often similar to each other, and can we harness this similarity to achieve speedup? Towards this end, we propose a phylogeny-like data structure, nc-tree, to compactly represent an ensemble of related biological networks; and a versatile framework, ProcessNets, that runs incremental algorithms over the nc-tree to speedup property computations on the input networks. Our asymptotic analysis of space and time complexity, coupled with empirical evaluations of diverse ensembles of simulated and real-world GTEx gene coexpression networks, demonstrate that ProcessNets leads to significant space and time advantages when computing different network science measures on ensembles of similar networks. For instance, we achieve a compression factor of 2.1x and 8.1x in storing two representative ensembles of GTEx Muscle Skeletal networks with 1000 nodes/genes, and a time speedup of 2.86x to 3.97x over the fastest baseline for computing degree centrality in ensembles of networks with at least 14,500 nodes derived from subsampled GTEx datasets of ten tissues. These results are promising and encourage application of ProcessNets to analyze biological network ensembles derived from rapidly accumulating consortium/biobank-based datasets.

In network biology, analyzing pathways is critical to study mechanisms or sequences of biological events. The pathway co-membership approach aims to link proteins belonging to the same pathway inside a network that is a clique. Although this method is great for projecting pathway knowledge onto proteins, it remains unsuitable for some downstream applications relying on graph distances like random walks with restart. To address this limitation, we introduce a weighted pathway co-membership method that incorporates biological distances between proteins. We achieve this by developing an abstraction framework for Reactome pathways, which captures both their hierarchical organization and sequential relationships. By assigning weights to these abstractions, we define a scoring system for proteins within the pathway co-membership clique. We illustrate our approach using a Reactome pathway ("Signaling by TGF-beta receptor complex", R-HSA-170834), and explain how we generate the weighted co-membership cliques of the 29 Reactome top-level pathways. Next, we construct the entire weighted pathway co membership network of the Reactome knowledge base. Validation and evaluation of our method is performed using cancer hallmark gene sets and a random walk with restart experiment, comparing our approach to state-of-the-art methods. We show that our approach outperforms state-of-the-art methods for gene prioritization on cancer hallmarks.

Next-generation sequencing data is compositional in nature and ignoring this fact can significantly impact downstream data analyses. This is especially true in co-occurrence (correlation) networks because compositional data is prone to many false negative correlations. The centered log-ration (clr) transformation is a commonly used compositional data analysis (CoDA) approach that converts compositionally constrained data into real number space. The clr transformed values can subsequently be used in machine learning and statistical analyses such as Spearman or Pearson correlations used to generate co-occurrence networks. Here we show that the necessary zero substitution performed prior to clr transformation can create numerous multicollinearity artefacts. We also show that a newer CoDA approach called the robust centered log-ratio (rclr) eliminates these artefacts in correlation-based network analysis.

We have developed a method to integrate transcriptomic coexpression networks across diverse experimental conditions within a single species. Our framework expands previous pannetwork approaches into a hypergraph based pannetwork. It first identifies coexpressed gene clusters within each individual dataset and extracts frequently co-expressed gene sets across multiple datasets using a frequent itemset mining algorithm. This process yields a hypergraph where each hyperedge is assigned a frequency (universality, U). We then extracted a subnetwork comprising high U hyperedges as the core network and formed modules within it. We applied our method to 106 Escherichia coli transcriptome datasets from the GEO database. The modularity of the core network peaked at a universality cutoff of 15, which was subsequently used to define it. Approximately 70% of the resulting core modules correspond to operons, and conversely, approximately 70% of all operons are covered by these core modules. We visualized the core network via an inter-modular network and analyzed core module dataset relationships using a modularity profile matrix. Based on these analyses, we successfully visualized the dynamic reorganization of the coexpression network in response to environmental changes in bacteria.

AO_SCPLOWBSTRACTC_SCPLOWShapley-value-based explanations, notably SHAP (SHapley Additive exPlanations), have gained prominence as a principled game-theoretic framework for local explanations and global feature importance. While exact Shapley value computation is exponential in feature count, TreeExplainer exploits the recursive structure of decision trees to achieve polynomial-time computation for tree-based ensembles. In many scientific applications, however, features are naturally organized into a priori groups reflecting domain knowledge, requiring explanations both across and within groups. The Owen value extends the Shapley value through a two-stage allocation rule that incorporates group structure while preserving fairness properties; yet, efficient algorithms for its computation remain limited. In this paper, we propose exact and Monte Carlo algorithms for computing Owen values in tree-based ensembles by combining hierarchy-guided group aggregation with tree-aware dynamic programming. The exact algorithm computes Owen values without sampling under the path-dependent characteristic function, which approximates the conditional expectation, whereas the Monte Carlo algorithm provides a scalable approximation that is unbiased for any prespecified sampling budget and converges almost surely as the sampling budget increases. We also provide global importance measures and visualization tools for structured, multi-resolution explanations. The proposed algorithms and tools are collectively referred to as TreeOwen. Through simulation experiments, we demonstrate the numerical accuracy and substantial computational gains of TreeOwen. We illustrate its practical utility using immunotherapy metagenomic data, showing how microbial genera (groups) and species (features) contribute to patient recovery.

Underlying the development of multicellular organisms is the process of cell differentiation, which is governed by the concerted and sequential change in gene expression. Various methods have been developed that employ scRNA-seq data to infer the position of a cell along a pseudo-temporal axis and identify relevant genes involved in the process. These trajectory inference methods typically rely on global transcriptomic changes and mathematical methods. However, overemphasis on large-scale transcriptomic changes may impair sensitivity to identify branching points and convergent trajectories, which are rather governed by small-scale transcriptional events. Motivated by this, we developed Moirai, a graph-based trajectory inference method that identifies gene expression patterns that change dynamically over a developmental continuum and leverages these to define a common pseudotime axis between all cells. In doing so, Moirai shifts the focus to individual gene dynamics, which enhances its ability to detect putative branching points that are masked by global transcriptomic similarities. We apply Moirai to four developmental datasets, where we demonstrate its ability to recover gene expression patterns of genes with a known involvement in the respective developmental process, motivating their use for defining a cells pseudotime. We furthermore show that Moirai can robustly infer gene expression patterns across different embedding approaches, highlighting the value of moving the focus of the inference process to the small-scale transcriptional dynamics.

Estimating differential gene expression is a common task in RNA-seq. Current methods mostly rely on normalization to reduce variance and increase accuracy. These methods are widely used and provide invaluable information about transcriptomic changes between biological conditions. However, widely used normalization methods are known to distort estimates of transcript differential expression when a majority of genes are upregulated or downregulated, or when the total RNA content per cell changes. Despite the presence of global expression shifts in a number of contexts, few methods exist that can provide accurate normalization and estimate linear models under this context without spike-in controls. Here, we present \textbf{GLORB} (\textbf{G}\textbf{L}\textbf{O}bal-shift \textbf{R}obust \textbf{B}ayesian model), a method for estimating generalized linear models under global upregulation. We develop two models that are able to recover differentially expressed genes and linear model coefficients with lower distortion of results. We show that our model's method of accounting for library size variance is consistent with DESeq2's median of ratios, and edgeR's trimmed mean of M-values under conditions when a minority of genes are upregulated and outperforms them under circumstances when most genes are either increased or decreased between groups. Finally, because our method does not rely on calculating geometric means for each gene it is able to work in datasets with much higher sparsity.

Single-cell sequencing provides a rich source of data that can be used to model and estimate the evolution of somatic cells within individuals. Many common models for modeling DNA evolution in phylogenetics are based on continuous-time Markov processes that operate on aligned sites. For haploid DNA, such as found in mitochondria, there are many well established models for these processes but since nuclear DNA is diploid, containing copies of both the maternal and paternal genome, it is natural to develop continuous-time Markov processes adapted to diploid DNA. One possible model is that both the maternal and paternal evolve independently according to haploid models so that, in principle, DNA comes with a phase. i.e. a label that tells if the site was derived from the maternal or paternal sequence. However, this information may not be available in practice and only an unphased diploid genotype is available. That is, an unphased model groups both of the possible maternal-paternal labels into a single state. This grouping process, called ``lumping'', is used by other methods as well and results in a coarser process. Care must be taken so that the coarser process can still be modeled as a continuous-time Markov chain; the lumped process of the original phased Markov chain may or may not be Markovian, depending on restrictions imposed on the generator matrix. Here, we explore the mathematical formulations of phased and unphased Markov processes, as well as a collection of other diploid DNA models from the literature.

Discovering gene regulatory networks (GRNs) from single-cell RNA sequencing (scRNA-seq) data is critical for understanding cellular function. Still, existing methods are limited by strong theoretical assumptions or high computational complexity. We introduce a multiple testing framework for GRN construction using Chatterjee's rank correlation coefficient, a nonparametric measure of dependence. Our approach overcomes the limitations of traditional methods while offering a transparent, scalable, and computationally efficient alternative to recent black-box machine learning models. Crucially, to address the non-independence of cellular observations inherent to scRNA-seq, we develop a data-driven algorithm for estimating robust testing cutoffs. Furthermore, we exploit the asymmetric nature of Chatterjee's correlation to propose a new test for active regulation, enabling the construction of biologically meaningful and directionally informed GRNs. We demonstrate that our method matches or outperforms state-of-the-art approaches in recovering true gene-gene dependencies and directed regulatory interactions from both simulated and real datasets, particularly for complex, non-linear dependencies, providing a powerful tool for dissecting complex GRNs.

Learning from imbalanced data remains a major challenge for graph neural networks (GNNs), as minority nodes are not only rare but also structurally marginalized within the graph. We address this issue with CLARA, a hierarchical learning framework that decomposes node classification into two stages: a coarse subgraph-level classifier that selects regions likely to contain minority instances, followed by a fine-grained node-level predictor within these regions. This design improves sensitivity while maintaining scalability, filtering out irrelevant areas and focusing learning on topologically meaningful neighborhoods. Experiments on benchmark graph datasets demonstrate substantial gains over established imbalance-handling methods, with CLARA reaching an F1-score of 88.3\%. The same strategy achieves significant improvements in protein-ligand binding site prediction, underscoring its broad and consistent effectiveness across both biological and general graph learning tasks.

Single cell RNA sequencing (scRNAseq) has revolutionized cellular biology over the past decade. At its most basic, a scRNAseq dataset is a cell by gene matrix with entries giving the transcription count of a particular gene in a particular cell. A common workflow applies log-normalization and z-scoring to the count matrix and then embeds the cells in a lower dimensional space using principal component analysis (PCA). A large number of computational tools use the PCA embedding as a starting point, but statistical models for the distribution of cells in the PCA space are lacking. As a result, many tools depend on ad-hoc assumptions that likely lead to signal loss. In this work, we address this gap using the metacell model, a statistical model of count matrices introduced by previous authors. We show that application of the workflow to count matrices generated under a metacell model yields PCA embeddings distributed as a Gaussian mixture model (GMM). To our knowledge, this is the first result connecting a statistical model of counts to the distribution of cells in PCA space. Our derivation depends on methods from random matrix theory (RMT) and we extend and establish connections between RMT and scRNAseq count matrices. We then use the predicted GMM to test the metacell model on 11 datasets. The predictions are accurate for counts generated by the metacell model but underestimate the variance of real cells for most of our datasets. We show that this error largely arises from gene correlations within metacells, which the metacell model does not capture. Some of our datasets are fit relatively well by the metacell model, suggesting that the fit of the metacell model to true datasets depends on sufficiently dense sampling. Finally, we show that the strength of noise in the PCA space varies significantly across metacells, suggesting that statistical modeling of PCA embeddings can improve existing downstream workflows. Our results provide a framework for deriving and testing statistical models of PCA embeddings from statistical models of counts and identify a limitation of the metacell description for scRNAseq datasets.

Accurate enzyme annotation remains a major bottleneck in translating rapidly growing protein sequence data into biological knowledge. Enzyme Commission (EC) prediction is particularly challenging because enzyme functions are organized hierarchically, annotations are often imbalanced across classes, and sequence similarity alone may be insufficient to resolve functional differences. To address these challenges, we developed ESM-ECForest, a two-stage framework that combines protein embeddings generated by the pretrained language model ESM-2 (Evolutionary Scale Modeling 2) with Random Forest classifiers. The first stage distinguishes enzymes from non-enzymes, whereas the second assigns one or more EC numbers to proteins predicted to be enzymatic. On an external benchmark comprising 25,778 protein sequences, ESM-ECForest achieved the highest weighted F1 score among the evaluated methods at all four EC levels, decreasing from 0.94 at Level 1 to 0.90 at Level 4. The largest relative improvements were observed for lyases (EC 4), ligases (EC 6), and translocases (EC 7), although EC 6 and EC 7 remained the most difficult classes internally. Visualization of the ESM-2 embedding space using Uniform Manifold Approximation and Projection (UMAP) revealed clustering patterns consistent with enzyme functional relationships, indicating that biologically relevant information is retained in the pretrained representations prior to supervised classification. These results support the use of pretrained protein language model embeddings as an effective foundation for enzyme annotation. By combining large-scale sequence representations with a lightweight supervised classifier, ESM-ECForest provides a scalable approach for EC prediction and may facilitate functional annotation of protein sequences derived from large genomic and metagenomic datasets.

Longitudinal non-Gaussian data are common in biomedical and microbiome studies, where responses may not be normally distributed. However, many existing dynamic graphical models are developed for Gaussian responses, while many non-Gaussian graphical models focus on cross-sectional dependence and do not explicitly estimate temporal transition effects. We propose a latent Gaussian dynamic graphical model for longitudinal non-Gaussian responses, including binary, count, and compositional outcomes. The proposed framework jointly estimates the contemporaneous conditional dependence network and the lag-one temporal dependence structure. We further develop a penalized Monte Carlo expectation maximization procedure for estimation. Simulation studies across multiple graph structures and outcome types show that the proposed method can recover both dependence structures, with more stable recovery for the lag-one temporal dependence structure and more structure-dependent recovery for the contemporaneous conditional dependence network. We apply the proposed method to longitudinal data from the Integrative Human Microbiome Project (HMP2) using a multinomial observation model. The results reveal persistent lag-one temporal dependence and contemporaneous conditional dependence among dominant microbial taxa, demonstrating the utility of the proposed framework for dynamic network analysis of longitudinal non-Gaussian biomedical data.

Below 30% pairwise sequence identity, alignment-based methods struggle to reliably distinguish true homologs from chance (Rost 1999), and enzyme function prediction degrades accordingly: on proteins in this regime, even advanced methods (CLEAN) achieves only 55.1% accuracy at full EC specificity on the CARE benchmark (Yang et al. 2024). To this end Protein Language Models (PLMs) have gained favor as alternatives. However, PLMs often encode only a subset of the biology (Heinzinger et al. 2024, Lin et al. 2023), whereas the understanding of enzyme function requires among other things a combination of sequence, structure and functional-context simultaneously (Ribeiro et al. 2023). In this work, we describe FuncSeek, a contrastive learning model which utilizes three diverse, complementary PLMs: ESM2 (to model evolutionary co-variation), ProstT5 (for bilingual sequence and structure embeddings), and ProteinBERT (for functional semantic similarities). Using SwissProt data, these 2816-D embeddings are labeled with Enzyme Commission numbers (EC) and are trained through a supervised contrastive head into a 256-D space. FuncSeek attains 64.6% nearest-neighbour EC4 accuracy on the CARE out-of-distribution benchmark set (ood30; proteins below 30% identity to training set), outperforming CLEAN (55.1%) and Diamond BLASTp (51.4%), and obtains 93.7% nearest-neighbour EC4 accuracy on the promiscuous, multi-functional enzymes benchmark (CLEAN, 69.4%). We also show that the learned representations transfer without retraining to the TrEMBL database, achieving 97.3% nearest-neighbour EC4 accuracy on a 8,031 BRENDA-validated enzyme set (Schomburg et al. 2004), never seen during training. Because only projected embeddings are stored in the target index, and function is inferred from an annotated reference set, we propose this paradigm for rapidly searching extremely large metagenomic databases, bypassing costly sequence alignment and annotation pipelines.

Model selection is a fundamental step in phylogenetic analysis that determines the best-fit model of sequence evolution for a given multiple sequence alignment. Popular model selection methods, such as ModelFinder, rely on statistical information criteria, such as the Bayesian Information Criterion (BIC) or the Akaike Information Criterion (AIC). However, these approaches are computationally expensive and the use of information criteria has been the subject of ongoing discussion. Recently, machine learning has emerged as a promising approach for phylogenetic model selection in both nucleotide and protein sequence analyses. ModelDetector is currently the only machine learning-based method for amino acid substitution model selection. However, because ModelDetector was trained on simulated data, it does not perform well on real datasets. Another limitation is that it does not support different rate heterogeneity across sites (RHAS) models. To overcome these limitations, we introduce ProtFinder, an efficient machine learning framework for protein model selection that predicts amino acid substitution models, RHAS models, and amino acid frequency models. To enable ProtFinder to work with real datasets, we employed a transfer learning strategy consisting of three stages: (1) initial training on large-scale simulated data, (2) joint training on both simulated and real data, and (3) final fine-tuning using real data only. Experimental results show that ProtFinder outperformed ModelDetector in amino acid substitution model selection. ProtFinder achieved comparable accuracy to the maximum likelihood method ModelFinder for substitution model selection on medium and large MSAs. It performs slightly better than ModelFinder in RHAS model selection and substantially outperforms it in amino acid frequency model determination. Notably, ProtFinder is up to 1,400 times faster than ModelFinder in terms of inference time, making it particularly suitable for medium and large datasets.

Single cell RNA-sequencing (scRNA-Seq) data are typically represented as cell-by-gene count matrices, which capture the expression of each gene as detected in the sampled cells; often a heterogeneous population of multiple different cell types or cell states. Almost all scRNA-Seq analysis workflows have a gene selection step prior to applying clustering algorithms which helps remove genes with low variability and hence reduce the high-dimensional gene space. A de-facto method for achieving selection of highly variable genes (HVG) uses dispersion and mean expression scores to evaluate the variability of each individual gene. However, methods based on direct mean-to-variance relationship for gene selection often suffer from susceptibility to variance instability and arbitrary determination of the optimal number of genes to use in downstream analysis tasks, additionally, they often prioritize genes with low abundance but high variance. Here, we propose an innovative method for selecting highly variable genes that is not based on mean to variance ratios: "Principal Genes (PG)" method; it utilizes the rotations (or loadings) from Principal Component Analysis (PCA) to calculate a novel variability score per gene that we name "Gene Principal Score (GPS)". GPS helps evaluate the genes based on their contribution in the PCA rotations and hence ranks the genes according to their variability from highest to lowest variable genes. For efficient implementation we utilize Augmented Implicitly Restarted Lanczos Bidiagonalization methods to efficiently obtain Principal Components (PCs) associated with the largest variance. Genes with the highest GPS score, i.e. Principal Genes, can then be used for downstream analysis tasks, especially the clustering step. To test the performance of our highly variable gene identification method, we use several validation strategies, including clustering of labeled single cell RNA-Seq data (i.e. data with known ground truth cell type labels). Furthermore, we measure the performance of our method against dispersion-based highly variable gene (HVG) selection approaches. We use several validation metrics, including sensitivity and adjusted rand index scores for clustering based on genes selected using our method against genes selected using HVG; and our validation datasets include six real labeled single cell RNA-Seq datasets. Our findings show that our new method, Principal Genes, is comparable and often favorable in performance in selecting highly variable genes and achieves ultra-fast gene selection from PCA results.

Retrieving BGCs related to those of a known producer can be regarded as a representation-learning objective. We hypothesize that ESM-2 sequence-derived representations of BGCs can improve retrieval beyond the Pfam-domain content metric. Our toolkit is the following: group-disjoint train, validation, and test assignments, validation-frozen model selection, [fi]ve seeds, and family-level paired inference. Of 6,953 atlas BGCs from 182 deduplicated Streptomyces griseus genome accessions, 5,325 silver-labeled BGCs are split into 98 training, 21 validation, and 21 test reference groups. Of the test reference groups, 16 are eligible for retrieval diagnostics. Pfam Jaccard scored Recall@50 of 0.8788, while Pfam-augmented BGC-SetNet scored 0.8472. The combination of ESM and Pfam-augmented BGC-SetNet scored 0.8769. A weighted Pfam Jaccard obtained a slightly higher score of 0.8789, which has a negligible difference compared to unweighted Pfam accard. Our results do not support the claim that sequence-derived representations can recover alternative biosynthetic pathways on this benchmark. Instead, explicit Pfam remains the major signal for this objective. Our results de[fi]ne the curation and pathway-level validation processes that are necessary for a more robust biological test.

BackgroundPhylogenetic tree reconstruction from single cell data based on copy-number alterations (CNAs) is an important problem in cancer genomics. Methods have been developed to address this problem by computing pairwise distances between copy-number profiles and employing a tree reconstruction algorithm. Despite the tight interplay between distance estimation and tree reconstruction, these two steps are often treated as separate problems, with the choice of the reconstruction algorithm receiving little attention. Most methods rely on classical Neighbor Joining (NJ), an algorithm designed for unrooted phylogenies that does not account for the fixed diploid root inherent to copy-number evolution. ResultsWe identify the Deepest Least Common Ancestor NJ (DLCA-NJ), not previously applied in this context, as the appropriate algorithm for phylogenies from copy-number data. By leveraging the known diploid root, it consistently outperforms standard NJ on simulated benchmarks across all evaluated metrics, with the most pronounced improvement in root placement accuracy. Building on these findings, we introduce Sparse Rooted Neighbor Joining (SRNJ), a scalable adaptation of DLCA-NJ. SRNJ significantly reduces running time while trading off only a minor loss in accuracy. We provide theoretical and empirical evidence of robustness to mutation rate using both synthetic and real biological datasets. ConclusionsRooted NJ variants offer a principled way to exploit the known diploid root when reconstructing phylogenies from copy-number data, and SRNJ extends this advantage to datasets whose size places the full distance matrix out of reach. The gains are clearest where distances are reliable, as on simulated data, while on real data accuracy appears to be constrained by distance estimation rather than by the reconstruction algorithm, leaving room for improvement as callers advance.

Gene perturbation experiments coupled with transcriptomic profiling are crucial for uncovering causal gene-gene relationships, yet it remains cost-prohibitive to systematically explore perturbation responses across diverse biological conditions. As a result, in silico prediction of perturbation response has emerged as an important strategy for guiding cost-effective experimental design. Although recent methods have begun to address cross-condition perturbation prediction, it remains under-characterized across scenarios defined by whether the perturbation has been observed during training under other biological conditions. Here, we study cross-condition prediction under both seen- and unseen-perturbation scenarios. We introduce TranScouter, a lightweight encoderdecoder framework that represents perturbed genes using LLM-derived embeddings of their text summaries and represents biological conditions using transcriptomic profiles of control cells from the target condition. Across evaluated benchmarks, TranScouter performs competitively in both scenarios. We further use empirical analyses to characterize how condition-space coverage and perturbation-effect transferability shape crosscondition performance.